# The following script was adopted from Thiele, Jan C., Kurth, Winfried and Grimm, Volker (2014) 
# 'Facilitating Parameter Estimation and Sensitivity Analysis of Agent-Based Models: 
# A Cookbook Using NetLogo and 'R'' Journal of Artificial Societies and Social Simulation 17 (3) 11 
# <http://jasss.soc.surrey.ac.uk/17/3/11.html>.

library(RNetLogo)
library(sensitivity)
library(lhs)

set.seed(1395878419)

# NetLogo path; NetLogo model path; simulation function path
nl.path     <- "C:/Program Files (x86)/NetLogo 5.2.0"
model.path  <- "C:/Users/wozni_000/Desktop/sobol/LMA2.0R.nlogo"
simfun.path <- "C:/Users/wozni_000/Desktop/sobol/sim_fun1_all.R"

# min. and max. values for the each of parameters      
input.values <- list(
  'search-bonus'  = list(quantile.function="qunif", min = 0, max = 5),
  'ltu-search-bonus' = list(quantile.function="qunif", min = 0 , max = 5),
  'agency-adds'    = list(quantile.function="qunif", min = 0,  max = 0.40),
  'ltu-agency-adds'          = list(quantile.function="qunif", min = 0, max = 0.4),
  'tau-regular'    = list(quantile.function="qunif", min = 0,  max = 0.3),
  'tau-ltu'      = list(quantile.function="qunif", min = 0, max = 0.3)
#  'benefits'      = list(quantile.function="qunif", min = 0.3, max = 1.1),
#  'minimum-wage'      = list(quantile.function="qunif", min = 1, max = 1.6)
)


# names of output values
output.names <- c("umean.criterion","theta.criterion","ltu.criterion", "wages.criterion", "ltu.wages.criterion", 
                  "duration.criterion", "ltu.duration.criterion")

# the number of repititions for the eFAST method
no.repititions <- 150

# how many repetitions for each input factor set should be run 

no.repeated.sim <- 10

plot.on.screen <- TRUE

# should R report the progress
trace.progress = TRUE

# initialize NetLogo
nl.LaborMarketACE <- "nl.LaborMarketACE"
NLStart(nl.path, gui=TRUE, nl.obj=nl.LaborMarketACE)
NLLoadModel(model.path,nl.obj=nl.LaborMarketACE)

# load the code of the simulation function
source(file=simfun.path)

# get names of input factors
input.names <- names(input.values)

# get names of quantile functions fpr the input factors
q.functions <- sapply(seq(1,length(input.values)), function(i) {
  input.values[[i]]$quantile.function})

# generate a list of arguments for the quantile functions
q.args <- lapply(input.values, function(i) {
  i$quantile.function <- NULL; return(i)})

# create instance of fast99 class
f99 <- fast99(model = NULL, factors = input.names, n = no.repititions, q = q.functions, q.arg = q.args)

# get required number of simulation iterations
iter.length <- nrow(f99$X)

already.processed <- 0

# simulation results for input factor sets (as matrix)
sim.results.f99 <- apply(f99$X, 1, function(x) {simulate(param.set=x,                            
                                                         no.repeated.sim=no.repeated.sim, 
                                                         nl.obj=nl.LaborMarketACE, trace.progress=trace.progress,
                                                         parameter.names=input.names, iter.length=iter.length,
                                                         function.name="eFAST")}
)

# iterate over evalulation criteria, print and plot results
plot.fast99 <- function(x, ylim = c(0, 1), ...) {
  if (! is.null(x$y)) {
    S <- rbind(x$D1 / x$V, 1 - x$Dt / x$V - x$D1 / x$V)
    colnames(S) <- colnames(x$X)
    bar.col <- c("red","blue")
    barplot(S, ylim = ylim, cex.axis=1.2, col = bar.col, names.arg=c(1,2,3,4,5,6), cex.names=1.5)
    legend("bottomleft", x.intersp=0.1, text.width=c(0, 2.8, 0.4, 0.6), c("first-order effects", "interactions efects"), 
           cex=1.2, fill = bar.col, bty="n", horiz=TRUE, xpd=TRUE, inset=c(0, -0.12))
    
  }
}

titles = c("unemployment density", "tightness fluctuation", "ltu density", "wages", "ltu wages", "duration", "ltu duration")


par(mfrow=c(4,2), mar=c(4,4,2,1), mgp=c(2,0.8,0), las=1)
for (i in (1:7))
{
  tell(f99, sim.results.f99[i,])
  
  if (plot.on.screen == TRUE) {
    par(ask=TRUE)
  }
  
  # plot results 
  print(f99)
  plot(f99, ylim=c(0,1))
  title(main=titles[i], cex.main=1.5)
  title(xlab="parameters", cex.lab=1.5, line=3.1)
  title(ylab="effects", cex.lab=1.5, line=2.5)
}     
